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The equation 

d t u = ud 2 xx u - (c - l)(d x u) 2 

is known in literature as a qualitative mathematical model of some biological phe- 
nomena. Here this equation is derived as a model of the groundwater flow in 
a water absorbing fissurized porous rock, therefore we refer to this equation as 
a filtration-absorption equation. A family of self-similar solutions to this equa- 
tion is constructed. Numerical investigation of the evolution of non-self-similar 
solutions to the Cauchy problems having compactly supported initial conditions is 
performed. Numerical experiments indicate that the self-similar solutions obtained 
represent intermediate asymptotics of a wider class of solutions when the influence 
of details of the initial conditions disappears but the solution is still far from the 
ultimate state: identical zero. An open problem caused by the nonuniqueness of 
the solution of the Cauchy problem is discussed. 
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1 A derivation of the filtration-absorption equation 



A permeable rock layer on a horizontal impermeable bed is considered. It is well known 
that the gently sloping groundwater flows in a purely porous medium with infiltration are 
described by the Boussinesq equation (see refs. (1-3)): 

mdth = n{hd 2 xx h + (d x h) 2 ) +q. [1] 

Here h is the groundwater level, t is the time, x is the horizontal space coordinate along 
the impermeable bed, q is the intensity of the groundwater inflow or outflow. Furthermore, 
k = pgk/ p is a joint property of the pair-rock/fluid, assumed here to be constant, k is the 
rock permeability, m the rock porosity, g the gravity acceleration, p and p are fluid density 
and dynamic viscosity. 

Assume now that when the fluid level is decreasing, some part of the fluid remains being 
absorbed by rock, e.g. due to capillary imbibition to the micropores. Then for a fixed fluid 
particle the rate of fluid absorption q will be proportional to the individual time derivative 
of the fluid level dh/dt, so that 

q = am ^ = am{d t h + u p d x h) . [2] 

Here a is also a joint rock/fluid property which we also assume here to be constant, and u p 
is the actual fluid velocity. 

According to the Darcy law the filtration velocity (fluid flux per unit area) is equal to 

The actual fluid velocity is different (see (1)) from the filtration velocity Uf. For a purely 
porous medium u p = Uf/m. It is very important that if the rock is fissurized, i.e. contains 
a connected network of cracks, then 

u P = ^ [4] 
mi 

where m x ( "fissure porosity") is much less than block porosity m: fluid is contained in pores 
but moves through cracks which are much wider than the pores but occupy much less space. 
Therefore 

q = am{d t h-^(d x h) 2 ) . [5] 
pm\ ' 

Substitution of [5] to the water balance equation [1] gives the equation 

m{\ - a)d t h = ^(hd 2 xx h + (1 - a—) (d x h) 2 ) , 
p m\ 
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which can be reduced to the form 



dth = m K ^ hd ** h - ( c " l )^ h f) . [6] 
m(l — a) 

where 

c = a — . 7 
mi 

Here there is a specially important point to be mentioned: If m\ = m (non-fissurized purely 
porous medium), then c = a and obviously c is always less than one because the absorption 
cannot exceed the available amount of fluid. However, if the fluid is fissurized, mi can be 
substantially less than m, and c = am /mi can be substantially larger than one. Replacing 
x hy x/ a/ K,/m(l — a) and leaving previous notation x for transformed space coordinate, 
we reduce the basic equation to the canonic form 

d t h = hd 2 xx h-(c-l)(d x h) 2 , [8] 

which will be investigated further. 

Note that the proposed model has some common features with the Mirzadzhan-Zadeh 
model (4) of filtration of the gas-condensate mixture, but is not identical to this model which 
leads to a different basic equation. 



2 Self-similar solutions 

We look for self-similar solutions with shrinking support and with finite time to collapse 
(total annihilation): For sufficiently large c (in fact, see later, for c > 3/2) 

Here xq is the point where the solution is collapsing at t — to- We assume the structure 
of the groundwater "dome" is symmetric, so that F(£), £ = (x — xo)/B(t — t) M is an even 
function. Here A, B and A are constants, but as we will see later, of a different nature. We 
will determine the function in the interval < £ < 1, so that = at £ > 1. The 
quantity 

x f = B(t -ty [10] 

is the contracting half-width of the groundwater dome. As is well known (5) for the de- 
generate parabolic differential equations of the type under consideration, the support of the 
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solution remains compact if it is compact initially. Furthermore, the function F can be 
normalized arbitrarily, so we can assume A = B 2 fi. Substitute [9] into the basic equation 
[8]. Bearing in mind that the coefficients of the resulting equation for F(£) cannot contain 
the time t explicitly, we obtain A = 2/i — 1, and the equation for assumes the form: 

We turn now to the boundary conditions. The first one is the condition of symmetry 

F'(0) = . [12] 

The second boundary condition follows from the continuity of the groundwater level at 
the free boundary x = Xf. 

F(l) = . [13] 

The last boundary condition follows from the continuity of the groundwater flux at the free 
boundary x — Xf. The solution close to free boundary can be considered as quasi-steady: 
h = h(Q , ( = x — Xf. We obtain from [8]: 

Tr dh , d 2 h . f dh\ 2 

where V = dxf/dt. 

For the solutions of equation [11] having the quantity d(F 2 )/dt; equal to zero at £ = 1 — 0, 
which is needed to have continuous flux at the free boundary, the quantity hd 2 h/d( 2 tends 
to zero at x = Xf — 0, so that at x — Xf 

Bearing in mind that for the self-similar solution [9], V = dxf/dt = —Bfi(to — t)^ 1 , and 
(dh/d()^ =0 = (d x h) x=Xf , we obtain the relation 

Bn{t Q - t)^ 1 = -Bfji{t - ty-\c - l)F'(l) 

from which the third boundary condition follows, to be satisfied by the function F 

*"(!) = _J_. [16] 
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So, a nonlinear eigenvalue problem is obtained: We have to find for the second order equation 
[11] in the interval [0, 1] the solution satisfying three boundary conditions [12], [13], [16], and 
the eigenvalue //. 

We consider the special case when the function F(£) has a maximum at £ = 0. This 
allows one to search the solution in the form of an expansion 

oo 

F(0=a(l-£ 2 ) + ^a n (l-£ 2 r , 

n=2 

so that the terms of the sum (except the first one) do not contribute to all three boundary 
conditions of the eigenvalue problem. The result is unexpectedly simple: 

1 c - 1 

a = — , a = , a n = (n > 2) [17] 

2(c- 1) 2c -3 ' v - ; 

so that, if c > 3/2, the solution to the nonlinear eigenvalue problem is obtained in the form 

1 'i-a, »=^\ [is] 



2(c- l) v * ' ' ^ 2c -3 
and the self-similar solution under consideration is represented by the relation 



h = whw) B2{t °- t)T ° 



_ (_X - Xq) 



2(c-l) 



[19] 



3 Investigation of the self-similar solutions 

In spite of its very simple form, the solution [19] is a typical self-similar solution of the 
second kind (see (6)). The exponent fj, — (c — l)/(2c — 3) cannot be obtained using some 
conservation laws, but only by solving a nonlinear eigenvalue problem, and the constants B 
and to ) as well as xo , are obtained by matching the self-similar solution with the solution to 
the Cauchy problem at the non-self-similar stage. To a certain extent this problem is similar 
to the problem of the evolution of a turbulent burst (see (6), section 10.2.4). The solution 
[19] has essentially different behavior in various intervals of the values of the absorption 
coefficient c: 

< c < 1; 1< c < § ; | < c . [20] 

The form [19] is appropriate for the last interval 3/2 < c where the collapse time t is finite. 
It is instructive to investigate the limiting behavior of the solution [19] at c — > 3/2 from 
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above. Putting c = 3/2 + e (e > is a small parameter), we obtain ji — I fie + |, so that 



and 



(t - *r = tt 



x f (t) = Btg 



+ \ 4 £ " r 2 
1-- 



{to - 1) 



j_ 

+ 2 £ 
''O 



1 - - 



1 — 

to 



— +± 

t \ 4e + 2 



[21] 



Therefore if at £ — > the quantity 4et tends to a certain constant G, and the quantity 
— +i 

B 2 t§ to another constant which we denote by C 2 9, the solution [19] tends to a finite 
limit: 



h = C 2 e~ 2t / @ 



2 1 



1 - 



C 2Q g-2t/G 



In the interval 1 < c < 3/2 the exponent /i becomes negative, and it is convenient to 
replace /i by —(/,, and t by — 1 . Solution [19] may be represented in a different form 



h 



-B 2 {t + ty — 



i - 



{X - XqY 



2(c-l) 
3-2c 



2(3 -2c)' [_ 5 2 (t + t)^ 

so that h{xo,t) = h mSuX {t) and x/ decay with time according to the power laws 

1 



[23] 



hraax(t) 



-B\t + t)- , x f (t)=B(t + t)- . 



[24] 



2(3 -2c)' 

The time of collapse is infinite and to becomes simply an additive constant. In the limit 
c — > 3/2 from below, c = 3/2 — e, e > 0, we obtain /i = l/4e — 1/2, and 



hmax(t) 



^bh;^ +1 

Ae 



tr 



1 



Bt 4 ^ 2 



1 + 



t 



tr 



tends 



Assuming again that at e — > the quantity 4et tends to a constant 6 , and £?t 4 
to another constant C\/©, we obtain the same limiting formula [22]. 

In the interval < c < 1 (weak absorption) the compact support extends, not contracts, 
although slower than in the case of "porous medium equation" c = 0. In this special case 
c = solution [19] is reduced to a known self-similar solution fo the first kind ((7), (8); see 
also (9), (6)). The degenerate special case c = 1 was considered previously; the papers by 
J. R. King (10) and P. Rosenau (11) should be mentioned specially. It is instructive to 
compare the results obtained above with those obtained in the paper by B. Meerson et al 
(12). 
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4 Nonuniqueness of solutions of the Cauchy problem 



We consider solutions of equation [18] with initial condition 



h(x, 0) = h (x) 



for x G R . 



[25] 



where ho(x) is a continuous function which is positive in an interval (xl(0),xr(0)) and 
which vanishes elsewhere. Let XL(t) and x R (t) be two continuous functions for t > such 
that Xi{t) is nondecreasing, x R (t) is nonincreasing, and XL(t) < xji(t) for t > 0. It is 
known that if c > 1 for any such pair x^(t) and xn(t) there exists a solution h(x,t) of the 
Cauchy problem [8], [25] such that h(x,t) is positive if ^(t) < x < xn(t), t > and h(x,t) 
vanishes elsewhere. For the proof, the definition of solution, and further references we refer 



Of special interest is the choice of steady interfaces: x L (t) = x L (0) and x R (t) = x R (0) 
for all t > 0. The corresponding solution is larger than any other solution, and in (14) 
a numerical scheme has been introduced which leads to this unique solution. From the 
modeling point of view it is interesting to observe that this solution can be obtained from 
the following limiting procedure: replace ho(x) by ho(x) +e (e > 0), solve problem [8], [25], 
and let e — > 0. 

Angenent (15) has constructed a solution of [8], [25] if h (x) has nonzero slope at xl(0) 
and xr(0) (for technical reasons the construction is local in time), which is unique in the 
class of solutions which can be expanded in a Taylor series of suffcieintly high degree near 
the interfaces: 



a similar expression holds at the right interface x = x R (t). Here the uniqueness not only 
refers to h, but also to the interfaces. We observe that the self-similar solution [19] belongs 
to this class of solutions. 

In section 5 we shall construct a numerical scheme which yields solutions converging to 
the self-similar solution for large times, and it is natural to ask if these solutions belong 
to the class introduced by Angenent. In particular this scheme yields solutions which are 
different from the ones obtained by the scheme in (14). We conjecture that the solutions 
which we construct in the following section are physically relevant, but undoubtedly future 
research is needed to provide definite answers to the nonuniquness question. 



to (13). 



N 



h(x, t) = J2 c feW(^ - x dt)t + o((x - x L (t)) N ) 



as x — > xiit) ; 
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5 Numerical experiment 



The goal of the numerical experiment was to indicate that the self-similar solution obtained 
above attracts the solutions to non-self-similar Cauchy problems having the initial condition 
of compact support, generally speaking a non-symmetric one: 



h(x, 0) = h (x) , x L (0) < x < x R (0) , 

and h(x, 0) = outside the interval xl(0) < x < x R (0). The basic equation 
transformed to a form convenient for numerical calculations 



can be 



Xp(t) -Xf(t)£ 1 

dth ~ */(t) d * h + m? 



[hd^h- (c- l)(d^hf 



where xq = dxo/dt, if = dxf/dt, and 



x - x (t) 



x f (t) 



x R (t) - X L (t) 



x (t) 



x R (t) +x L (t) 



so that the interval of new space variable £ becomes fixed: — 1 < £ < 1, whereas the 
solution h(x,t) is different from zero in the time dependent interval x L {t) < x < x R (t). 
Again, assuming naturally the quasi-steadiness of the level distribution in the vicinities of 
free boundaries, we obtain the conditions 



[26] 



at t > 0, and the basic equation takes the form: 
1 



d t h = 



x){t) 



[27] 



with the initial condition 



Me,o) = VO» iei<i; MC,o) = o, |£| > 1 . [28] 

Two numerical schemes were used in our computations performed by finite-difference 
approximations: (i) a forward-in-time, centered-in-space explicit approximation, and (ii) 
a forward-in-time, centered-in-space implicit approximation. For the most part numerical 
calculations have been run with the time step At = 1CT 5 for the explicit scheme and At = 
10~ 4 for the implicit one. The number iV of subintervals of length A£, iV = 2/A£ was equal 
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to 202 for both schemes. The results obtained by using these numerical approximations 
coincided with good accuracy. The absorption coefficient c was always equal to 1.75. 

The first initial condition was taken as a "smoothed block": a homogeneous water level 
distribution smoothly going to zero at the edges. The results of the computation are pre- 
sented in Figure 1 in the form of the distribution of the scaled level: level h(x, t) divided by 
maximum level at each time h max (t). It is seen that the curves corresponding to different 
times collapse to the parabola, corresponding to the self-similar solution [18]. The time of 
collapse t and the constant B were determined in the following way (Figure 2): According 
to the intermediate asymptotics [19] at small t — t, 

/wx(t) ^(0) • [ J 

so that the quantity a^(t)//i max (t) should be a linear function of time, and the intersection 
of its graph with the time axis (Figure 2, a) gives the value of t . We obtain from [10] a linear 
relation between \nxf and ln(t — t), i.e. in the coordinates — \nxf, — ln(t — t) a straight 
line with the slope /i. It gives us the value B (Figure 2,b) and an additional possibility of 
checking the asymptotics. Naturally to an d B depend on the initial condition. For our case 
we found to = 0.345, B = 3.73 and ji = 1.499, which agrees well with the analytic value 
H = 1.5. 

The next computation was performed for a nonsymmetric initial condition: 

-Ax 2 + Ax , < x < \ 



h(x, 0) 



- gX + gX + g , 2 ^ *^ 2 . 



Figures 3, a and 3,b demonstrate the behavior of the numerical solution for different 
times. It is clearly seen that the solution becomes symmetric and tends to the self-similar 
asymptotics [10]. The values of to — 1-138, B = 0.63 and /i = 1.5027 have been calculated 
as before; the calculated and analytic values of ji agree with high precision. 

For comparison we have taken the solution [19] for a certain t as an initial condition, and 
computed the solution to the partial differential equation further using the same algorithm. 
The results are presented in Figure 4, a for different times. Being plotted in scaled coordinates 
(Figure 4,b) they collapse to a single curve, giving us an additional check of the numerical 
procedure. 

6 Conclusions 

We presented a new derivation of the filtration-absorption equation based on a model of 
groundwater flow with partial absorption. It is shown that for a sufficiently large absorption 
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constant the time of collapse is finite. A family of self-similar solutions to this equation 
is obtained. Numerical experiments indicate that these self-similar solutions obtained are 
self-similar intermediate asymptotics for the solutions to the Cauchy problems having the 
initial conditions data with compact support, but due to the nonuniqueness of the solution 
of the Cauchy problem future research is needed to provide more definite conclusions. 
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Figure Captions 



Figure 1. The numerical solution to the Cauchy problem for c = 1.75 with the initial 
condition of a "smoothed block" type for different times in the scaled coordinates. The 
solution is collapsing to the self-similar asymptotics. 

Figure 2. 

(a) The determination of the parameter t - 

(b) The determination of the parameter B. 

Figure 3. 

(a) Evolution of the nonsymmetric initial distribution to a symmetric self-similar 
asymptotics. 

(b) The same evolution presented in scaled coordinates. 
Figure 4. The numerical solution preserves the self-similarity. 
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